TV-body Simulations for f(R) Gravity using a Self-adaptive Particle- Mesh Code 
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We perform high resolution N-body simulations for f(R) gravity based on a self-adaptive particle- 
mesh code MLAPM. The Chameleon mechanism that recovers General Relativity on small scales is 
fully taken into account by self-consistently solving the non- linear equation for the scalar field. We 
independently confirm the previous simulation results, including the matter power spectrum, halo 
mass function and density profiles, obtained by Oyaizu et al. (Phys.Rev.D 78, 123524, 2008) and 
Schmidt et al. (Phys.Rev.D 79, 083518, 2009), and extend the resolution up to k ~ 20 h/Mpc for 
the measurement of the matter power spectrum. Based on our simulation results, we discuss how 
the Chameleon mechanism affects the clustering of dark matter and halos on full non-linear scales. 
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I. INTRODUCTION 

The late-time acceleration of the Universe is the most 
challenging problem in cosmology. Within the framework 
of general relativity (GR), the acceleration originates 
from dark energy with negative pressure. The simplest 
candidate for dark energy is the cosmological constant. 
However, in order to explain the current acceleration of 
the Universe, the required value of the cosmological con- 
stant must be incredibly small. Alternatively, there could 
be no dark energy, but a large scale modification of GR 
may account for the late-time acceleration of the Uni- 
verse. Recently considerable efforts have been made to 
construct models for modified gravity as an alternative 
to dark energy, and to distinguish them from dark energy 
models by observations (see [1-5] for reviews). 

Although fully consistent models have not been con- 
structed yet, some indications of the nature of the modi- 
fied gravity models have been obtained. In general, there 
are three regimes of gravity in modified gravity models 
[2, 6]. On the largest scales, gravity must be modified sig- 
nificantly in order to explain the late time acceleration 
without introducing dark energy. On the smallest scales, 
the theory must approach GR because there exist strin- 
gent constraints on the deviation from GR at solar system 
scales. On intermediate scales between the cosmological 
horizon scales and the solar system scales, there can be 
still a deviation from GR. In fact, it is a very common fea- 
ture in modified gravity models that gravity can deviate 
from GR significantly on large scales. This is due to the 
fact that, once we modify GR, there arises a new scalar 
degree of freedom in gravity. This scalar mode changes 
gravity even below the length scale where the modifica- 
tion of gravity becomes significant that causes the cosmic 
acceleration. 
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In order to recover GR on small scales, it is re- 
quired to suppress the interaction arising from this scalar 
mode. One of the well known mechanisms is the so-called 
Chameleon mechanism [7]. A typical example is f(R) 
gravity where the Einstein-Hilbert action is replaced by 
an arbitrary function of Ricci curvature (see [8, 9] for a re- 
view). This model is equivalent to the Brans-Dicke (BD) 
theory with a non-trivial potential. The BD scalar me- 
diates an additional gravitational interaction, which en- 
hances the gravitational force below the Compton wave- 
length of the BD scalar. If the mass of the BD scalar 
becomes larger in a dense environment like in the solar 
system, the Compton wavelength becomes short and GR 
can be recovered. This is known as the Chameleon mech- 
anism [7]. In the context of f(R) gravity, by tuning the 
function /, it is possible to make the Compton wave- 
length of the BD scalar short at solar system scales and 
screen the BD scalar interaction [10-13]. 

This mechanism affects the non-linear clustering of 
dark matter. We naively expect that the power-spectrum 
of dark matter perturbations approaches the one in the 
ACDM model with the same expansion history of the 
Universe because the modification of gravity disappears 
on small scales. Then the difference between a modified 
gravity model and a dark energy model with the same ex- 
pansion history becomes smaller on smaller scales. This 
recovery of GR has important implications for weak lens- 
ing measurements because the strongest signals in weak 
lensing measurements come from non-linear scales. One 
must carefully take into account this effect when con- 
straining the model in order not to over-estimate the de- 
viation from GR. 

Using a concrete example in f(R) gravity, Refs. [14-16] 
performed N-body simulations to confirm this expecta- 
tion. They showed that, due to the Chameleon effect, the 
deviation of the non-linear power spectrum from GR is 
suppressed on small scales. It was also shown that the fit- 
ting formula developed in GR such as Halofit [17] failed 
to describe this recovery of GR and it overestimated 
the deviation from GR. Refs. [14-16] used a multi-grid 
technique to solve the non-linear equation for the scalar 
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field. To make predictions for weak lensing, it is nec- 
essary to model the power spectrum on smaller scales. 
Ref. [18] derived a fitting formula for the non-linear power 
spectrum in the so-called Post Parametrised Friedmann 
(PPF) framework [6] based on perturbation theory and 
N-body simulations. Ref. [19] calculated weak lensing sig- 
nals for the future experiments using this fitting formula. 
It was shown that the constraints on the model crucially 
depend on the modeling of the non-linear power spec- 
trum. 

The resolution of N-body simulations performed in 
Refs. [14-16] is limited to k < a few h Mpc -1 for the 
power spectrum mainly due to that fact that they used 
a fixed grid size. In this paper, we exploit a technique 
to recursively refine a grid to solve the scalar field equa- 
tion and the Poisson equation aiming to probe the power 
spectrum down to k ~ 20 h Mpc -1 . We modify the pub- 
licly available N-body simulation code MLAPM and add a 
scalar field solver. We also study properties of halos in 
our simulations. The Chameleon mechanism depends on 
local densities thus its effect depends on the mass of ha- 
los. It is important to quantify the effectiveness of the 
Chameleon mechanism to maximise our ability to distin- 
guish between models. 

This paper is organised as follows. In section II, we de- 
scribe our f(R) gravity models and present basic equa- 
tions to be solved. In section III, we explain the details of 
our N-body simulations and how we analyse the simula- 
tion data to obtain the power spectrum and halo proper- 
ties. In section IV, we present the matter power spectrum 
and properties of halos. We first use a conservative cut-off 
for the matter power spectrum to check the consistency 
of our results with linear predictions on large scales and 
previous results in Refs. [14-16]. Then using the high 
resolution simulations, we present the power spectrum 
on small scales and develop the fitting formula that cap- 
tures the effect of the Chameleon effects based on the 
PPF formalism. We then study the properties of halos by 
measuring the mass functions and halo density profiles. 
We also measure the profiles of the scalar field inside ha- 
los to quantify the effects of the Chameleon mechanism. 
Section V is devoted to conclusions. 



II. f(R) GRAVITY AND CHAMELEON 

In this section we briefly summarise the main ingredi- 
ents of the f(R) gravity theory and its properties, which 
are essential for the simulations and discussions below. 
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in which G is the Newtonian constant and Cm is the 
Lagrangian density for matter fields (radiation, baryons 
and dark matter, which in this paper is assumed to be 
cold) . Taking a variation of the above action with respect 
to metric yields the modified Einstein equations 

G a p + X a p = 87rGT a p. (2) 

Here the extra term X a p denotes the modification of GR, 

X a (3 = f R R a (3 - ( 77 - J 9afi D a Dpf R , (3) 



where f R = denotes an extra scalar degree of free- 

dom, dubbed scalaron, D a denotes the covariant deriva- 
tive with respect to metric and □ is the Laplacian oper- 
ator. Taking the trace of the modified Einstein equation 
Eq. (2) yields 

Uf R = ^f = \{R - f R R + 2/ - SnGp) = 0, (4) 

which governs the time evolution and spatial configura- 
tion of f R . Given the function form of f(R), a complete 
set of equations governing the dynamics of the model is 
obtained. 

We shall work here in an almost Friedman- Robertson- 
Walker (FRW) universe with the line element given as, 

ds 2 = a 2 (rj) {[1 + 2$(£,7?)]d77 2 - [1 - 2tf(f, rj)dx 2 }} , 

(5) 

in which r] is the conformal time, a(rj) is the scale fac- 
tor, is the gravitational potential and \£ is the spatial 
curvature perturbation. Subtracting off the background 
part of this equation, we obtain a dynamical equation for 
the perturbation of the scalaron, 



V 2 <% = - — [6R(f R ) + SirGSpu], 



(6) 



where Sf R = f R (R)-f R (R),SR = R-R, Sp M = Pm~Pm, 
V is the three dimensional gradient operator and we work 
in the quasistatic limit, meaning that we assume the spa- 
tial variation of the scalaron is much larger than the time 
variation, so that we could neglect the time directives of 
f R and approximate DSf R as V 2 5f R 1 . Here the quan- 
tities with a bar denotes those evaluated in the cosmo- 
logical background. Notice that Sf R , 5R and SpM are not 



A. The f(R) Model 

A possible generalisation of GR is to add a term, which 
is a function of the scalar curvature R, to the Einstein- 
Hilbert action S (see [8, 9] for a review and references 



1 The validity of this assumption has been verified by Ref. [14], 
namely, they found that the amplitude of the 2-norm of the spa- 
tial derivative is larger than that of the time derivative by a 
factor of 10 5 at all redshifts. We also checked the validity of the 
quasistatic assumption with our simulations, and found a good 
consistency with Ref. [14]. 
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necessarily small, and we call them perturbations just for 
convenience. 

We can also obtain the counterpart of Poisson equation 
for the gravitational potential by adding up the 00 and 
ii components of the modified Einstein equation Eq. (2) 
(see Appendix A of [20] for some useful expressions using 
the above metric convention) 



WttG 



a 2 Sp M + —6R(f R ), 
3 6 



(7) 



in which we have neglected terms such as 4> and as 
we are working in the quasistatic limit, and used Eq. (6) 
to eliminate V 2 S/r. 

Given the density field, a functional form of f(R) and 
the knowledge of the background evolution, Eqs. (6) and 
(7) are closed for the scalaron and the gravitational po- 
tential. These are thus the starting point for our TV-body 
simulations. 

Note that the remaining modified Einstein equations 
give a relation between two metric perturbations $ and 



V 2 (^ - $) = V 2 



(8) 



where we assumed Jr <C 1 in the background. Then com- 
bining Eqs. (6), (7) and (8), we find that the relationship 
between the lensing potential ^ + <£> and the matter over- 
density is the same in both GR and f(R) scenarios, which 



is, 



V 2 0£ + $) = 87rGa 2 5p : 



The Chameleon mechanism 



(9) 



It is well known that f(R) gravity models can incor- 
porate the so-called Chameleon mechanism [7, 21], which 
are vital for their cosmological and physical viability 
[10, 22]. The Chameleon mechanism was first proposed 
in the context of coupled scalar field theories, and used 
to give the scalar field an environment-dependent effec- 
tive mass so that is very heavy in dense regions 
and thus the scalar-mediated fifth force is suppressed lo- 
cally so as to avoid conflicts with experiments and ob- 
servations. Because the f(R) gravity is equivalent to a 
coupled scalar field theory in the Einstein frame with the 
functional form of f(R) related to the potential of the 
scalar field, a similar mechanism is needed to ensure that 
the fifth force in f(R) gravity is well within the limits set 
by experiments. 

To see how the Chameleon mechanism works in f(R) 
gravity models, let us consider their difference from GR. 
For Eq. (6), 5/r vanishes identically in GR, and we have 
SR = — SttGSpm which, substituted into Eq. (7), recovers 
the usual Poisson equation for GR, 



V 2 $ = 4irGa 2 5p. 



(10) 



Therefore, to make Chameleon mechanism work, we have 
to ensure that in Eq. (6) V 2 5Jr = V 2 Jr is close to zero 2 , 
meaning that /r must depend on R very weakly (i.e., be 
nearly constant) in high curvature region. Furthermore, 
if fn deviates significantly from zero, then Eq. (3) means 
that the theory will become one with a modified Newto- 
nian constant, which is not of our interest. A necessary 
condition for the Chameleon mechanism to work here is 
thus | | <C 1 at least in the high curvature region 3 . 

The functional form of f(R) must be carefully designed 
so that it can give rise to the late time cosmic acceleration 
of the universe without any conflict with the solar system 
tests by virtue of the Chameleon mechanism. One such 
model was studied in [22], where f(R) oc (—R) n with 
n <C 1 (n = corresponds to a cosmological constant). 
A more interesting model was proposed in [10], namely, 



f(R) 



-m 



2 ci(-#/ra 2 ) n 
c 2 (-R/m 2 ) n + 1' 



(11) 



where m 2 = 87rGpM,o/3 = HqQm, note the minus sign 
in front of R is due to our different sign convention from 
[10]. To evade the solar system test, the absolute value 
of the scalaron today, |/ro|j should be sufficiently small. 
However the constraint is fairly weak (|/ho| < 0-1) [10] 
and it is satisfied in our models. In this scenario, the 
scalaron always sits at the minimum of the effective po- 
tential defined in Eq (4), thus, 



R « 8ttGpm ~ 2/ = 3m 2 ( a~ 3 + - — V 

V 3c 2 / 



(12) 



To match the ACDM background evolution, we need to 
have 



Cl _ 
C2 



(13) 



Plugging in the numbers w 0.76, ~ 0.24 into 
Eq (12), we have — R « 41m 2 ^> m 2 , which can be used 
to simplify the expression for the scalaron, 



ci n(-i?/m 2 ) n " 1 
"4l(-R/m 2 ) n + 1] 2 



/ m 2 y 



(14) 



Therefore we can see that the two free model parameters 
are n and ci/c 2 , and the latter is related to the value of 
the scalaron today via 



ci 



; ( 1 + 4 fT 



71+ 1 



#0, 



(15) 



2 This means that Sfn is only sensitive to the local matter density 
distribution, which is merely a reflection of the fact that the extra 
scalar degree of freedom 5/r has a heavy mass (in dense regions) 
and cannot propagate away from the source (and as such the fifth 
force is severely suppressed). 

3 There is also requirement on the sign of fji to avoid instability 
in the perturbation growth [22, 23]. 



We shall concentrate on the models with n = 1 and 
\Iro\ = 10" 4 , 1(T 5 , 1(T 6 in this paper, as in [14]. 

In the linear regime where SpM/pM <C 1, we can lin- 
earised the scalaron equation by linearsing Eq. (12) with 
respect to a cosmological background 

Sf R = -(n+l)f B0 (5° J +1S -§, (16) 

where R is the background curvature and the subscript 
implies that the quantity is evaluated today. Then the 
scalaron equation (6) becomes 

V 2 % = aV%-^«. ( 17 ) 

where 

In Fourier space, the solution can be easily obtained as 

1 87rGa 2 5p M nm 
5fR= 3k* + a*P' (19) 

Then it is easy to show that we recover GR above the 
Compton wavelength a/k > X c as S/r is suppressed by 
the mass term compared with the gravitational potential, 
\^Ir\ ^ 3>- This is because the scalaron cannot propa- 
gate beyond the Compton wavelength. Below the comp- 
ton wavelength, 5/r = —87rG5py[/3 and we get 

V 2 $=^a 2 ^ M , (20) 

which implies that the gravitational constant is enhanced 
by 4/3. This leads to the scale-dependent enhancement 
of the linear growth rate. This linearisation around the 
cosmological background fails when Sf R becomes larger 
compared with the background field Sf R ^> \/r\. This 
happens in the dense region because S/r is driven by 
the matter perturbations SpM- In these dense regions, 
the curvature is large and Eq. (16) ensures that /r is 
suppressed and GR is recovered, realising the Chameleon 
mechanism. In Fig. 1, we plot the time evolution of the 
Compton wavelength and the background /r field with 
n = 1 and \f R0 \ = 1(T 4 , 1(T 5 , 1(T 6 , which is useful to 
understand the recovery of GR in our simulations. 

III. THE iV-BODY SIMULATIONS 

In the model of [10], the fifth force due to the propaga- 
tion of S/r is suppressed by many orders of magnitudes 
than gravity only in very dense regions, while in less dense 
regions it has similar strength as gravity. Such a dramatic 
change of the amplitude of the fifth force across space in- 
dicates that Eq. (6) should be highly nonlinear, which 



100 




redshift z 



FIG. 1. The evolution of the Compton wavelength A c (upper 
panel) and the absolute value of the background scalar field 
as a function of redshift z. The f(R) models of |/ro| = 
10" 4 ,10" 5 ,10" 6 are illustrated using black solid, red dashed 
and green dash-dot curves respectively. 



the readers can convince themselves by appreciating the 
nonlinearity of Jr [cf. Eq. (14)]. Consequently treatments 
involving linearisation fails to make correct predictions, 
and TV-body simulations which explicitly solve the non- 
linear equation for S/r are needed. 

In a series of papers [14-16], Oyaizu et al. pioneered 
in this direction, by solving Eq. (6) on a regular mesh to 
compute the total force on particles. Their results show 
some interesting and cosmologically testable predictions 
of the f(R) gravity. However, the resolution of their sim- 
ulations was limited by their regular mesh. Here, we per- 
form similar simulations, but with an adaptive grid which 
self-refines in the high density regions. This technique has 
been applied previously in [24-26] and works well. 

As we are interested in how well the Chameleon mech- 
anism actually works, besides the full f(R) simulations, 
we also perform some non- Chameleon runs with the 
Chameleon effect artificially suppressed (see below). For 
clarity we shall refer to these two classes of simulations 
respectively as full f(R) and non-Chameleon simulations. 



A. Outline of the Simulation Algorithm 

We have modified the publicly available TV-body sim- 
ulation code MLAPM [27], which is a self-adaptive particle 
mesh code, for our TV-body simulations. It has two sets 
of meshes: the first mesh includes a series of increasingly 
refined regular meshes covering the whole cubic simula- 
tion box, with respectively 4, 8, 16, • • • , cells on each 
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side, where Nd is the size of the domain grid, which is the 
most refined of these regular meshes. This set of meshes 
are needed to solve the Poisson equation using multigrid 
method or fast Fourier transform (for the latter only the 
domain grid is necessary). When the particle density in 
a cell exceeds a predefined threshold, the cell is further 
refined into eight equally sized cubic cells; the refinement 
is done on a cell-by-cell basis and the resulted refinement 
could have arbitrary shape which matches the true equal- 
density contours of the matter distribution. This second 
set of meshes are used to solve the Poisson equation using 
a linear Gauss-Seidel relaxation. 

Our core change to the MLAPM code is the addition of 
the routines which solve Eq. (6), which is similar to that 
in [25, 26]. These include: 

1. We have added a solver for S/r, based on Eq. (6), 
which uses a nonlinear Gauss-Seidel relaxation iter- 
ation and the same convergence criterion as the de- 
fault Poisson solver of MLAPM. It adopts a V-cycle in- 
stead of the default self-adaptive scheme in arrang- 
ing the Gauss-Seidel iterations, because the latter is 
found to be plagued by the problem of oscillations 
between adjacent multigrid levels, and the conver- 
gence is too slow. 

2. The value of S/r solved from the above step is used 
to complete the computation of the source term for 
the Poisson equation Eq. (7), which is solved us- 
ing a fast Fourier transform on the regular domain 
grids and Gauss-Seidel relaxation on the irregular 
refinements. 

3. The value of solved from the above step is used 
to compute the total force on the particles by a fi- 
nite difference, and then the particles are displaced 
and accelerated as in usual TV-body simulations, ac- 
cording to 



(21) 
(22) 



where a subscript c denotes code unit. More details about 
the code could be found in [26, 27]. 



B. iV-body Equations in Code Units 

To implement Eqs. (6, 7) into our numerical code, we 
have to rewrite them using code units, which are given 

by 

X p 

X c = — , Pc = TT j-O ^0 = tflQ, 





- E£ 


dt c 




dp c 


1 


dt c 


a 



B 1 



tt^, Pc = — , V c = BV, (23) 
Pm 



where B is the size of the simulation box and Hq = 
100ft- km/s/Mpc. We also need to choose a code unit 
for S/r (or but this is done in different ways for 
the Chameleon and non- Chameleon simulations for con- 
venience, and here we shall discuss them separately. 



1. The Full f(R) Simulations 

In our full f(R) simulations, the quantity SR in Eq. (6) 
is written in its exact form SR = R — R with 



R 



-m 



£1 J_ 

^ 2 £ 



R = —3m 



2 ' a" 3 



(24) 
(25) 



where Q m and Q\ are respectively the present fractional 
energy densities for matter and dark energy, and m 2 = 
QmHq. With the aid of these, we can rewrite Eq. (6) 
using code units as 



ac 



(BH ) 



R 



i of ci n " 

-tt M (p c - 1) + -n M a \-n— — 



"4 Sr 



l + 4^a 3 



(26) 



Note that we will drop the subscript of V c for simplicity 
from now on. 

In the simulations, |/#| spans a range of many orders 
of magnitude, from ~ 10 -13 in the very dense regions to 
~ 10 -5 in the low density regions. Obviously, the above 
equation is very sensitive to numerical errors in the trial 
solution to /r, making it difficult to solve Jr accurately. 
Furthermore, Eq. (24) indicates that /r>0 would lead 
to unphysical (imaginary) R, and if this happens by ac- 
cident in the simulations due to the numerical errors 
(which is very likely because \/r\ <C 1) then the com- 
putation cannot carry on further. To solve these prob- 
lems, we introduce a new variable u = log(— /r) (which 
is slightly different from the choice of [14]). Then, typi- 
cally u G (—30,-10), which does not vary much across 
the simulation box and Jr = — e u < no matter which 
value u takes. 

Replacing /r with u, we can rewrite Eq. (26) as (after 
some rearrangement), 



vV = n M p c + 40v^ 



(BHoY 
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n M ^(n^) e -^+r. (27) 



Similarly, in terms of u and using the code units, we 
can rewrite the Poisson equation as, 

3 

V 2 $ c = 2n M p c - -n M + 2^ A a 3 



(HoB) 2 



l^u {^-^j + e" 7 *^ 3 - ( 28 ) 



6 



6 



Eqs. (27, 28), after discretisation (Appendix A), can 
be implemented in the modified MLAPM code straightfor- 
wardly. 



2. The N on- Chameleon Simulations 

As mentioned above, because the Chameleon effect is 
one of the most important features of our f(R) model, 
we also would like to see what happens if we artificially 
suppress it. Unlike [15, 16], in our simulations we achieve 
this by explicitly linearising the equation for Jr with- 
out using the solution for S/r given in Eq. (19). Since 
the Chameleon effect originates from the nonlinearity of 
fn, the linearisation will largely weaken it, justifying the 
name non- Chameleon simulations. 

The linearisation is easily done by simply setting SR = 
~§j^3 which changes Eq. (6) to 



c 2 V 2 f R = c 2 V 2 



8ttG 



Sp M a 



1 OR 



Sf R a 2 . (29) 



3df R 

In this case, we find it more convenient to choose 5/r as 
our variable, then using the relation 



OR 



,2°2 



1 



a 



n+2 



c\ n(n + 1) [ V 
we can rearrange the above equation to obtain 

~2 



(30) 
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1 flMd 3 



2 V z 5f R 
-1) 



3 n(n + 1) c\ 



n+2 



Sf R . (31) 



Note that after the linearisation the mass of the scalaron 
depends only on the background matter density; but if 
\fRo\ is chosen to be small enough (such that c?>/ci ^> 1) 
the mass can be very heavy, and thus the fifth force still 
gets suppressed, even though the suppression is not as 
strong as in the full f(R) simulations. Furthermore, the 
suppression of the fifth force in this case is universal and 
independent of the local matter density, again justifying 
its name - non-chameleon model. 

Similarly, the modified Poisson equation appears to be 

V 2 3> 



2n M ( Pc -i) 

1 £1m& 3 < 



6 n(n + 1) c\ 



n+2 



Sf R . (32) 



Eqs. (31, 32), after discretisation (Appendix A), can 
be implemented in the modified MLAPM code directly. 



C. Code Tests 

As the major modification to the MLAPM code is for 
the S/r solver, we have to check carefully whether it 





Box size (Mpc/h) 
256 128 64 


N ■ 

1 v sim 


10 


10 


10 


N P 


256 3 


256 3 


256 3 


iVgrid 


128 


128 


128 


&Ar/2(n/ M P c ) 


0.79 


1.57 


3.14 


/c*(h/Mpc) 


5.5 


11.0 


22.0 


Refinement levels 


8 


9 


10 


Force resolution (Kpc/h) 


12 


23 


94 


Mass resolution (lO n M /h) 


13.3 


1.75 


0.21 



TABLE I. Technical indices for the simulations presented in 
this work. The quantity iV s i m denotes the number of reali- 
sations we run for each model, N p is the total number of 
particles, iVg r id is the number of domain grids one each side, 
and k N /2 and which are defined in Eq. (41), show the half 
Nyquist scale and the /c-cutoff we actually used in this work 
respectively. 



works correctly. For this purpose we have solved Eq. (27) 
around a point-like mass [14] and compared the result 
with the analytic solution; we find that the numerical and 
analytic solutions agree as well as in [14] (Fig. 2 there). 
To avoid repetition we shall not display the results here. 

Other, indirect, tests of our code include comparing the 
predicted matter power spectrum to the linear perturba- 
tion results to see whether they agree on large scales. We 
shall discuss these in turn as the paper unfolds. 



D. Simulation Details 

We choose to run our simulations in a cosmology that 
is consistent with the WMAP observations. Specifically, 
we parametrise the universe in the form of 



p = {c,q}. 



(33) 



The symbols C and Q denote the subsets of cosmological 
and f(R) parameters, respectively. 

C = {n h h 2 , Q c h 2 , Sl K , H , n s , a 8 }, (34) 

where ^t>^ 2 and fl c h 2 denote the physical energy density 
for baryons and cold dark matter respectively, Qk is the 
curvature, Ho is the Hubble constant today, n s stands 
for the index of the primordial power spectrum at a pivot 
scale of k = 0.05 Mpc -1 , and as measures the amplitude 
of the (linear) power spectrum on the scale of 8 Mpc/h. 

g = {n,f R0 }. (35) 

Here n and f R o specify a f(R) model as described in 
Sec. II. We set the cosmological parameters as C = 
{0.04181, 0.1056, 0, 73, 0.958, 0.8}, while for the f(R) pa- 
rameters, we fix n = 1 and simulate for three models with 
fm = -10" 4 ,-10" 5 ,-10" 6 . We also run the standard 
GR model (/r = 0) for the purpose of comparison. For 
each f(R) model, we also run the simulation without the 
Chameleon mechanism using the same initial condition 
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and background evolution. Since we will refer to these 
simulations frequently later in the text, we assign them 
abbreviations to lighten the notations, as listed below, 

• F4 = Full f(R) simulation with \f R0 \ = 1(T 4 ; 

• F5 = Full f(R) simulation with \f R0 \ = 1(T 5 ; 

• F6 = Full f(R) simulation with \f R0 \ = 10~ 6 ; 

• N4 = non-Chameleon simulation with \f R o\ = 
10- 4 ; 



• N5 = non-Chameleon simulation with \f R o\ = 
10- 5 ; 

• N6 = non-Chameleon simulation with \f R o\ = 
10- 6 ; 

• GR = Simulation for a ACDM model in GR. 

For each model listed above, we run 10 simulations 
with 10 different initial conditions (ICs), which we call 
10 'realisations', to reduce the sample variance. To ex- 
tend the range of scales, we run these simulations in 3 
difference boxes with size B = 256, 128 and 64 Mpc/h. 
The technical details are summarised in Table I. We gen- 
erate the ICs using Grafic, an IC generator included 
in the COSMICS package [28], at redshift z = 49 in GR, 
and then evolve the system using our modified version 
of MLAPM for f(R), as well as using the default MLAPM for 
GR simulations. 



E. Data Analysis 

Data analysis is of great significance for the simulation. 
In this section, we will detail our pipeline to obtain the 
snapshots, matter power spectra, mass function and the 
density and scalar field profiles out from the raw simula- 
tion output. 



1. Snapshots 

Visualisations of the simulation result is helpful to un- 
derstand the physics in an intuitive way. For this pur- 
pose, we will show the 2-D snapshots for the distribution 
of overdensity S = 5p/p, the perturbation of the scalar 
field Sf R and gravitational potential <£>. For the visualisa- 
tion, we output the data for ln(l + 5), —2Sf R and $ on a 
400 x 400 x 400 grid from our simulation, and project the 
3-D volume onto a 2-D plane to make image snapshots. 
The resolution we use here is much lower than that we 
use for calculation in the code, but it is sufficient for the 
purpose of visualisation. 



2. Binning, Average and Variance 

We will present our simulation results (power spectra, 
mass function and profiles) in terms of data bins along 
with error bars. Suppose our observable is called O which 
is a function of x, and we have iV raw raw samples mea- 
sured from simulation for O in some range of x that we 
are interested in, then we take an average of O over 10 
realisations first, make logarithmic bins in x, and then 
calculate the mean value 0{ and error bar cr(O^) for the 
ith bin via, 

Oi = Ys°ij/ N ^ °\O t ) = Y J (0 %3 -6 t ) 2 /N l , (36) 

where Ni is the total number of raw samples falling into 
the ith bin. 



3. Matter Power Spectrum 

The matter power spectrum P(k) is an important sta- 
tistical measure of the matter clustering on different 
scales in Fourier space. We measure the matter power 
spectrum using a tool called P0WMES [29]. P0WMES esti- 
mates the Fourier modes of a particle distribution based 
on a Taylor expansion of the trigonometric functions, and 
it is able to accurately determine and correct for the bi- 
ases induced by the discreteness and by the truncation 
of the Taylor expansion. Also, the aliasing effect is safely 
negligible if the order of Taylor expansion N = 3, which 
we chose to use. One could even accurately measure P(k) 
up to the scale of k ~ 10247r/ (Box size) with the 'folding' 
operation, which we adopted in our analysis. 

We firstly measure the power spectrum for GR simu- 
lations using P0WMES, and use Pgr as the observable O. 
Then we utilise Eq. (36) to estimate the mean value and 
error bar for Pgr(&) for each k bin and compare it to the 
Halofit (Smith et al., Ref. [17]) prediction. Then we mea- 
sure the power spectra for full f(R) and non- Chameleon 
simulations, use the relatively difference 

AP/P GR =(P f{R) - P G r)/P G r, (37) 

as the observable O and similarly use Eq. (36) to get 
estimates of AP/Pgr. 



4- Halo Mass Function 

We use the MHF tool [30] (the Halo-finder for MLAPM) 
to resolve halos in our simulation. The halo-finding algo- 
rithm we use is similar to that in Ref. [16] - we assign the 
particles to the grids using a Triangular Shaped Cloud 
(TSC) interpolation, and count the particles within a 
sphere around the highest overdensity grid point. We in- 
crease the radius of the sphere until the overdensity S 
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reaches a threshold A. This process can be repeated un- 
til all the halos are identified, and we only keep those 
that are made of at least 800 particles for our analysis. 
In this paper we choose A to be the virial overdensity 
in ACDM model, A = A v i r . Thus the mass is defined 
by M = 47rA v i r r^ ir /3 where r v i r is the radius when the 
overdensity reaches the threshold 6 = A v i r . A v i r depends 
on redshifts; A v i r = 373.76 at z = and A v i r = 242.71 
at z = 1 in our GR models. We use the same defini- 
tion of the mass in f(R) simulations in order to make 
a fair comparison between f(R) gravity models and GR 
models. However, we should bear in mind that in f(R) 
gravity the virial overdensity is different from that in GR 
[16] and a care must be taken when we compare the mass 
function in our f(R) simulations to observations. 

In this work we use the definition of the halo mass func- 
tion as the halo number density per logarithmic interval 
in the virial mass M in GR, i.e., 

n ]nM ss (38) 

This definition is different from Ref. [16] where they use 
A = 300 to define the mass. Thus a direct comparison is 
not possible without rescaling the mass. 

Since the halos in f(R) and GR simulations have dif- 
ferent number and mass in general 4 , and we are in- 
terested in the relative difference of the mass function 
in f(R) model with respect to that in GR, we need to 
make sure that we are comparing the same quantity in 
different gravity models. Therefore, for the simulations 
with the same box size, we find the overlapping halo 
mass range for full f(R), non- Chameleon and GR sim- 
ulations for all the 10 realisations, and make logarith- 
mic bins in mass in this range. Then we count the num- 
ber of halos falling into each mass bin for f(R) and GR 
cases respectively, and calculate the relative difference 

AninM/rcinM fafnM ^mV^IhM for each bin for eV " 

ery realisation. Finally we average over 10 realisations to 
calculate the mean value and variance for each mass bin. 



5. Halo Profile 

We modified MHF to analyse the profile of the pertur- 
bation of the scalar field 5f Rl and the gravitational po- 
tential as well as the overdensity 6 as a function of 
the rescaled virilised halo radius r/r v i r . To see the rela- 
tive difference of the density profile with respect to that 
in GR, we have to identify the same halos in the full, 
non-Chameleon and GR simulations for each realisation. 
Since more halos are produced generally in f(R) models, 
we start from the GR halos in each realisation, and de- 
cide whether to include it for the analysis in the following 
steps, 



4 The number and mass of the halos are in general different even 
for the same gravity model for different realisations. 



1. For the first realisation in GR, record the position 
and mass for the first halo, also the number of par- 
ticles made up of this halo; 

2. For the same realisation for the 6 f(R) simulations 
[F4,F5,F6 and N4,N5,N6], search for the same halo 
as it in GR. If all the following conditions are sat- 
isfied, 

(a) The number of particles made up of this halo 
should be greater than 800; 

(b) In all the 6 f(R) simulations, a halo with the 
similar position (absolute distance to the cor- 
responding halo in GR smaller than 0.1% of 
the box size) and a similar mass (relative mass 
difference \AM/M - 1| < 1000%) is found; 

we assume it is the same halo as in GR simulations. 

3. Repeat the halo identification process so that all 
the halos with counterpart in f(R) are selected; 

4. For every selected halo, rescale r to r/r v i r and lin- 
early interpolate the overdensity 6 to obtain 5 for 
the same r/r v i r for different gravity models; 

5. Calculate the relative difference AS/5gr = ($f(R) — 
£gr)/^gr and feed A5/5gr into Eq. (36) to obtain 
the statistics for density profile. 

It is interesting to test the working efficiency of 
Chameleon mechanism inside the halos, and this can be 
realised by comparing the profiles of Sf R to that of <E>. If 
Chameleon does not work at all where the mass term of 
the scalar field can be ignored, i.e., SR = 0, from Eqs. 
(6) and (7), we have 

- 25 f R = (39) 

Therefore we define an efficiency parameter T for the 
Chameleon mechanism as, 

r = | - 25f R /$.\ (40) 

If r ~ 1 then Chameleon hibernates, and if T <C 1 then it 
means Chameleon is very active. For the selected halos, 
we calculate T for both full and non- Chameleon simula- 
tions and feed it into Eq. (36) to obtain the statistics for 
the profile of T of the Chameleon mechanism. 



IV. RESULTS 

In this section, we will present our simulation results 
including the snapshots, matter power spectrum, the halo 
mass function and the profiles of the halo density and the 
efficiency parameter T. 
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FIG. 2. The snapshots for density perturbation 8, scalaron perturbation SfR, and the gravitational potential taken from 
the full f(R) simulations (6 left panels) and the non-Chameleon simulations (6 right panels) for three values of |/ro|, namely, 
\f R0 \ = 1(T 4 , 10" 5 , 10" 6 from left to right. The box size of the simulation is 64 Mpc/h, and all the snapshots are taken at z = 0. 
Note that we use a same color scheme for all the plots of —2dfR and namely, — 2(5/#,<E> G [—2.0,0.1] x 10 -5 , so that one 
could directly compare them to see the working efficiency of the Chameleon mechanism (see text for details). We use another 
color scheme for all the plots of the density distribution, say, ln(l + 6) £ [3.0, 8.0]. 



A. Snapshots 

First, let us view the snapshots first to get some basic 
idea on the physics of our simulations. In Fig. 2, we show 
the snapshots for the density perturbation <$, the scalaron 
perturbation 5f R , and the gravitational potential $ taken 
from the full f(R) and non-Chameleon simulations for 
three values of \f R o\. These snapshots are taken from 
one realisation with the box size 64 Mpc/h at z = 0. 
The density plots show the familiar web-like structures, 
while the scalar field and potential have similar pattern, 
but are smoother. From the snapshots we can see that 
—26 fn ~ $ for the F4, F5, N4 and N5 simulations, and 
the Chameleon works very well for F6 cases (\Sf R \ <C |3>|) 
These are natural - the Chameleon mechanism works 
when Sf R > \f R \ and Sf R < $ ~ 10" 5 . Thus the F5 
simulations are the critical case where the Chameleon 
mechanism is about to fail today. Note that in the N6 
simulation we also find that \25f R \ is smaller than |$|, 
because here the scalaron also has a heavy mass and a 
short Compton wavelength (though it is independent of 
local densities). 

In Fig. 3, we show the same quantities for the F4 
case at higher redshifts. Interestingly, we can see that 
the Chameleon does work very well at z > 3, and it hi- 
bernates at z < 3. This is because the background field 
| f R \ is small at high redshifts (see Fig. 1) thus the condi- 



tion for the Chameleon mechanism, Sf R > \f R \, is easily 
satisfied. 



B. Matter Power Spectrum 

In this section, we present the result for matter power 
spectrum in f(R) theory, in comparison with that in GR. 



1. Tests for GR simulations 

To start with, we test the accuracy of our GR simula- 
tions by comparing our results to the Halofit prediction 
[17]. The result is shown in Fig. 4. As one can see, our 
B = 128 Mpc/h simulation is very consistent with the 
Halofit prediction up to k ~ 10 h/Mpc, namely, the rel- 
ative difference is smaller than 10% on all scales. The 
B = 256 Mpc/h simulation tends to underestimate the 
power at k > 0.4 h/Mpc, which is due to the fact that 
we do not have enough particles to resolve small scales in 
a large box. The B = 64 Mpc/h simulation shows much 
more power than that predicted by Halofit on scales k > 5 
h/Mpc. This does not necessarily mean that our simula- 
tion is not accurate on those scales, instead, the Halofit 
prediction has never been tested on such small scales, 
thus it might be inaccurate. 



10 



+ 



e 



z=5 



z=3 



z=l 



z=0 



i 



■ 




'.41 









FIG. 3. The snapshots for density perturbation 5, scalaron perturbation 5 fit, and the gravitational potential $ taken from the 
full f(R) simulations with \/ro\ = 10 -4 at four redshifts, namely, z = 5,3, 1,0 as illustrated in the figure. The simulation box 
size and color scheme are exactly the same as that in Fig. 2. 



2. Low-resolution tests for f(R) 

After confirming that GR simulations are consistent 
with the Halofit predictions, we switch to the f(R) sim- 
ulations. We did the full f(R), Non-Chameleon and 
GR simulations in three different boxes with the size 
B = 256,128,64 Mpc/h, and assemble them to extend 
the range of scale. Generally, the P(k) result is reliable 
up to the scale of fc*, which is defined as, 

K = N eS xk N/2: k N/2 = 7rNy s /(W), (41) 

where iV e ff is a factor determined by the adaptive nature 
of the code, and iV e ff = 1 for the non-adaptive simula- 
tions, i.e., the pure particle- mesh (PM) simulations. The 
quantity k N / 2 1S the half Nyquist scale, and N p and B 
are the number of particles and the box size respectively 
(see Table I for the values we use in this work). Since 
for the first step we aim to reproduce the results in Ref. 
[16] up to their resolution (half Nyquist scale), namely, 
k < 3 h/Mpc, for a cross-check, we use the same fc-cutoff 
criterion as Ref. [16], namely, here we set 7V e ff = 1. We 
will present our results on smaller scales, which is still 
reliable, in the next section. Given the spectra for dif- 
ferent boxes, we take a volume- weighted average on the 
overlapped scales. In Fig. 5, we show the matter power 
spectra for both the full f(R) simulations and for the 



non- Chameleon cases at redshift z = 0. We took a ratio 
with respect to the GR power spectra with the same ini- 
tial condition, and then averaged over ten different reali- 
sations to reduce the sample variance, as in Refs. [14-16]. 

As one can see, the growth is generally enhanced in 
f(R) gravity, and the peak value of the enhancement 
increases as \/ro\. This is natural bacause, for large 
|/i*o|j e -9-, \Iro\ = 10 -4 , the Chameleon does not work, 
therefore the effective Newton's constant is enhanced 
by 4/3 below the Compton wavelength, resulting in a 
large enhancement of AP/Pqr. For small |/#o|> e.g., 
|/i*o| = 10 -6 , Chameleon works very efficiently, which 
suppresses the structure growth back to that in GR lo- 
cally, yielding a tiny enhancement of AP/Pqr. The case 
for \/ro\ = 10 -5 is critical - it is something between these 
two extreme cases, in which Chameleon works, but not 
as efficient as the |/ho| = 10 -6 case. Therefore the differ- 
ence between full and the non- Chameleon simulations in 
this case is intermediate. 

In Fig. 5, we overplotted the linear prediction, and the 
Halofit for f(R) gravity. The Halofit result is obtained by 
applying the Halofit fitting formula naively on the linear 
prediction for f(R). Our simulations are consistent with 
linear predictions on large scales. The Halofit prediction 
overestimates the deviation from GR for smaller \/ro\. 
This is because the Halofit uses only the information of 
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FIG. 4. Matter power spectra for standard GR at redshift 
z — (upper panel). The simulation results from different 
boxes are shown with error bars, while the Halofit result is 
over-plotted (green solid curve) for comparison. In the lower 
panel, the relative difference with respect to the Halofit result 
is shown, and the dashed magenta line illustrates AP — to 
guide eyes. 

linear power spectrum and fails to capture the effect of 
the Chameleon mechanism, which suppresses the devia- 
tion from GR. Note also the Halofit is calibrated in GR 
simulations. Thus once the deviation from GR becomes 
larger on smaller scales, we cannot trust the validity of 
the Halofit prediction. Our result is in prefect agreement 
with Refs. [15, 16]. 

3. High-resolution results for f(R) 

The self-adaptive nature of our code allows us to go 
beyond the half Nyquist scale, in other words, 7V eff can 
be greater than 1. The realistic value of N e ^ is largely 
determined by the maximum number of refinement trig- 
gered in the whole simulation process, which is 8 ~ 10 for 
our simulations. To be conservative, we choose 7V e ff = 7 
in all cases, and show the high-resolution result for power 
spectra in Figs. 6. As in Fig. 5, we show the relative dif- 
ference of power spectra in f(R) with respect to that in 
GR at various redshift s. 

Let us see the higher resolution result from B = 64 
Mpc/h simulations shown in Fig. 6 at z = 3, 1 and 0. 
Looking at the non-Chameleon result first, one could 
have the following observations, 



FIG. 5. The relative difference of matter power spectra in 
f(R) models with respect to that in GR at redshift z = 0. 
From upside down, we show the ratio AP/Pqr for both full 
f(R) (square) and for the non-Chameleon simulations (trian- 
gle) with | /ro I = 10 -4 , 10 -5 , 10 -6 . For comparison, we over- 
plot the prediction from Smith et al (blue solid) and from 
the linear perturbation theory (red dashed). The data points 
with error bars are assembled from results of different box 
sizes based on the volume-weighted average method with a 
conservative k cut-off. See text for details. 

it goes down after reaching a peak at z = and 
z = 1 for N4 and N5 cases; 

2. For the same model, the peak moves towards larger 
k as one goes to higher redshift; 

3. The ratio monotonically increases with k at z = 3 
for N4 and N5 cases and at all times for N6 cases. 

One could use GR as an example to understand these 
observations. It is true that GR is different from our non- 
Chameleon f(R) models, but the physical origin of these 
features are similar. In Fig. 7, we show the relative differ- 
ence of matter power spectra in GR with as = 0.9 with 
respect to that with erg = 0.8 at redshifts z = 3, 1 and 0. 
We use CAMB 5 to generate the linear power spectra, and 
use Halofit to account for the nonlinear it ies. 

Interestingly, the pattern here is similar to what we 
saw just now - the ratio in GR shows a peak, and as 
z increases, the peak becomes more pronounced, and it 
shifts to larger k. This is due to the transition between 



1. The ratio increases as one goes to small scales, and 



5 Available at http://camb.info 
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FIG. 6. The relative difference of matter power spectra in f(R) models with respect to that in GR at redshift z — 3 (left 
panels), z — 1 (middle panel) and z — (right panels) taken from B — 64 h/Mpc simulations. From upside down, we 
show the ratio AP/Pgr for both full f(R) (black squares) and for the non-Chameleon simulations (green triangles) with 
I /ho | = 10 -4 , 10 -5 , 10 -6 . In each panel, we show the Compton wavelength Ac and the absolute value of the background scalar 
field \f R \ 



the 2-halo term to the 1-halo term in the halo model de- 
scription of the power spectrum (for a review see [31]). 
If (j 8 is higher, the non-linearity becomes important ear- 
lier and the transition from 1 to 2-halo term appears at 
smaller k. Since the 2-halo term has a larger amplitude, 
this creates a peak in the ratio. As the transition wave 
number between 1 and 2-halo terms moves to smaller 
k for smaller z, the peak is shifted to smaller k. The 
same thing happens in non- Chameleon simulations. In 
this case, the linear power has the /c-dependent enhance- 
ment. The transition from 1-halo to 2-halo term appears 
at high k at z = 3 in f(R) models but we do not see the 
peak as the ^-dependent enhancement of the linear power 
spectrum hides the enhancement due the 2-halo term. At 
z = 0, the transition appears at smaller k where the lin- 
ear enhancement is weak then we see the peak due to the 



transition from 1 to 2-halo term. 

Having understood the non- Chameleon simulations re- 
sults, now let us look at the full Chameleon simulations. 
At z = 3, we can see that the Chameleon mechanism 
works well even in the F4 case, and the deviation from 
GR is effectively suppressed. Once the Chameleon does 
not work well, e.g., at z = for F4, the full simula- 
tion result traces the non-Chameleon result as expected, 
with slightly lowered amplitude. Interestingly when the 
Chameleon mechanism has started to fail, e.g. at z = 1 
for F4 and at z = for F5, the power spectrum in 
full simulations has larger amplitudes than those in the 
non- Chameleon simulations on small scales. This im- 
plies that complicated dynamics is taking place when the 
Chameleon mechanism fails to work. We will come back 
to this issue later by studying the density profiles of ha- 
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FIG. 7. The relative difference of matter power spectra in GR 
with as = 0.9 with respect to that with as = 0.8 at redshifts 
z — 3 (blue dash-dot), z — 1 (red dash) and z — (black 
solid). Halofft is used to calculate for the nonlinearities. 



los. For the case where the Chameleon mechanism works 
very well, i.e., the F6 case, the full simulation result has 
much less amplitudes, implying that GR is successfully 
recovered. 



4. PPF fit 

Hu and Sawicki proposed a fitting formula (PPF fit [6]) 
for the matter power spectrum in modified gravity based 
on the assumption that the power spectrum resembles 
that in GR on small scales, which takes the form of, 



P(k,z) 



•Pnon-GR(fc, z) + C„i£ 2 (fc, z)P GR (k, z) 



Y, 2 (k,z) 



(42) 



where Pgr denotes the non-linear power spectrum in a 
ACDM model with the same expansion history as that 
in the modified gravity models, and P n on-GR means the 
nonlinear power spectrum in modified gravity without 
the mechanism that recovers GR on small scales, i.e. the 
Chameleon mechanism. Therefore in our case, P n on-GR 
is exactly the power spectra for the no n- Chameleon sim- 
ulation, i.e., Pnon-GR = ^no-cham- E measures the non- 
linearities and once the non-linearities are significant 
E ^> 1, the power spectrum approaches Pgr- Koyama 
et at. [18] did a successful PPF fit and recovered the 
full f(R) simulation up to k ~ 0.5 h/Mpc presented in 
[15, 16] using Eq. (42) with E given by 



E 2 (M) 



2tt 2 



P\in(k,z) 



1/3 



(43) 



Redshift 




z 1 




z = 


l f l 
\jRo\ 


1 n~ 4 


in - 5 in-6 

1U 1U 


1 n -4 


in - 5 in - 6 
1U 1U 


C n l 


0.17 


1.1 3.4 


0.09 


0.44 1.6 


a 


0.42 


2.5 1.3 


0.26 


0.76 0.39 


P 


-0.22 


-2.2 fixed to 


-0.01 


-0.85 fixed to 


7 


1.6 


0.15 fixed to 


2.04 


0.24 fixed to 



TABLE II. The best-fit PPF parameters for various f(R) 
models at redshifts z = 1 and z = 0. 



Since our simulation goes to much smaller scales than 
[15, 16], we study whether we can extend this fitting for- 
mula to smaller scales. To be conservative, we use the 
results from B = 128 Mpc/h simulations to extend the 
formula up to k = 10 h/Mpc (see Fig. 8). We generalise 
Eq. (43) by adding three more parameters a, j3 and 7, 



E 2 (M) 



k 3 

^Plin(M) 



(44) 



where Pn n is the linear power spectrum in f(R). 



Combining Eqs. (42) with (44), we fit AP/P GR up to 
k ~ 10 h/Mpc, and the best-fit parameters are listed in 
Table II. The best-fit power spectrum curves for f(R) 
models are overplotted with the simulation data in Fig. 
9. As one can see, the agreement for F6 model is within 
one percent level even if we fixed f3 and 7 to zero. This 
implies that, if the Chameleon works, we can use the 
PPF formulae by varying c n \ and a only. However, if the 
Chameleon fails, the power spectrum tends to go back 
to the non-Chameleon spectrum on small scales, so that 
we have to introduce additional parameters f3 and 7 to 
parametrise this behaviour. Once we include these pa- 
rameters, we can fit the power spectrum in F4 and F5 
models accurately. In the quasi non-linear regime k < 1 
h/Mpc, we can ignore the corrections described by (3 and 
our results are roughly consistent with [18]. The PPF 
fitting formula would be useful when we try to find fit- 
ting fomulae similar to Halofit in f(R) gravity models 
because the non-Chameleon simulations are far easier to 
run because of the linearity of the equation and the cos- 
mological parameter dependence of the PPF parameters 
was found to be weak at least in the quasi non-linear 
regime. 



C. Mass Function 

1. GR 

We test our mass function in GR firstly by comparing 
it to the theoretical predictions proposed by Sheth and 
Tormen [32], and by Tinker et al. [33]. The comparison 
is shown in Fig. 10. We stack our results from different 
boxes to extend the mass range, and we see an agreement 
with theoretical predictions. It seems that our simulation 
underestimates the number of small halos (halos in the 
first bin) because of the lack of particles, but this problem 
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FIG. 9. The relative difference of matter power spectra in full f(R) simulations with respect to P(k) in GR at redshift z — 1 
(left panel) and z — (right panel). The results for three f(R) models \fiio\ = 10 -4 , 10 -5 , 10 -6 are shown in black blocks, 
red circles and green triangles respectively. The blue solid curves show the PPF fit, and the horizontal dashed line illustrates 
AP = 0. 
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FIG. 10. Upper panel: the mass function for GR at redshift 
z = stacked using B = 64 Mpc/h (black blocks), 128 Mpc/h 
(red circle) and 256 Mpc/h (green triangle) overplotted with 
the theoretical predications from Sheth-Tormen (black dash) 
and from Tinker et al. (blue solid); Lower panel: relative dif- 
ference of the simulation result with respect to the Tinker et 
al prediction. 

can be alleviated to some extent when taking the ratio 
of the mass function of f(R) with respect to that in GR. 

2. f(R) 

The result for the halo mass function is presented in 
Fig. 11. As one can see, our result at z = is consistent 
with that in Ref. [16] (cf Fig. 2). Let us understand our 
result by looking at the no n- Chameleon result first. In 
this case, the strength of gravity is enhanced by 4/3 be- 
low the Compton wavelength A c , regardless of the value 
of \fRo\. Meanwhile, A c increases with \/ro\, making the 
enhancement on the mass function more pronounced for 
large \/ro\ models. This is what we see in Fig. 11. Also, 
at redshift z = 0, the enhancement for larger halos is 
more significant, while at redshift z = 1, smaller halos 
are more populated. This is because the stronger gravity 
in f(R) models create more small halos at high redshift 
and assemble them to make more large halos at low red- 
shift compared with ACDM models. 

For the full f(R) simulations at z = 0, the number 
of massive halos is generally reduced compared to the 
non-Chameleon model predictions, because the suppres- 
sion of the fifth force due to the Chameleon mechanism 
is stronger in high density regions, where massive halos 
most likely reside. For the case of \/ro\ = 10 -4 , in which 
Chameleon does not work, the full and non-Chameleon 
results agree well, while for the \/ro\ = 10~ 6 case, the 



number of large halos is almost the same as the GR pre- 
diction due to the Chameleon effect. At redshift z = 1, 
the suppression starts from lower mass bins simply be- 
cause Chameleon was working more efficiently at higher 
redshift. 



D. Halo Profile 

We show the halo profile for the overdensity and for T 
at redshift z = 1 and z = in Figs. 12 and 13 respec- 
tively. We show the results for the halos in three mass 
ranges: [10 12 ,10 13 ), [10 13 , 10 14 ) and [10 14 , lO 15 ]M /h, 
taken from the B = 64, 128 and 256 simulation boxes 
respectively. 

A quick observation of the result is that on average, the 
profiles of the overdensity for f(R) gravity are similar to 
that in GR, namely, the relative difference is smaller than 
20% in all cases. This is consistent with the analysis in 
Ref. [16]. The complex patterns of the relative difference 
of f(R) and GR halo density profiles can be understood 
qualitatively as follows. 

1. For the F6 simulation, the Chameleon effect is very 
efficient throughout the cosmic history and so the 
fifth force has always been strongly suppressed, so 
that the predicted halo density profile should be 
very close to the GR result, which is what we see 
from Figs. 12 and 13 (the only exception is the case 
of low mass halos at z = 0, which is because these 
halos have lower density and so the Chameleon ef- 
fect is weaker, especially at late times). 

2. For the N6 simulation, the fifth force is only weakly 
suppressed (especially at late times), and the cen- 
tral attractive potential towards the halo centre is 
stronger than in GR, which means that particles 
tend to fall towards the halo centres, producing 
higher density profiles than the latter. 

3. The N5 and N4 results follow the same pattern as 
N6, but are more complicated as in some cases the 
density profile can be even lower than in GR. This 
is perhaps because the central attractive potential 
is not the only thing affecting the halo density pro- 
file. As the fifth force is unsuppressed from early 
times, the particles are typically faster than they 
are in GR, and tend to escape the halo, flattening 
the density profile. Note that fifth force both deep- 
ens the central potential and speeds up the parti- 
cles, and the latter is an accumulative effect, which 
makes it difficult to see which effect is dominating. 

4. The F4 simulation agrees with N4 simulation very 
well at late time, because Chameleon effect is in- 
significant. But at early times it produces higher 
density profile than N4, except for the small halos 
(for which Chameleon effect is again unimportant), 
which seems to indicate that the particles in the F4 
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FIG. 11. The relative difference of the mass function for f(R) models with respect to that in GR at redshift z — 1 (left) and 
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halos are still relatively slow, because the fifth force 
has not become strong enough for long. 

5. The case of F5 is further complicated by the fact 
that \/ro\ — 10 -5 is a critical point for our f(R) 
models. For the small halos, the fifth force has just 
become unsuppressed at z = 1, and we see that 
their density profiles are higher than N5 results for 
the same reason as above, while at z = the F5 and 
N5 simulations agree very well. For the medium- 
sized halos, at z = 1 the fifth force is still strongly 
suppressed and their density profiles are lower than 
N5 results because the central attractive potential 
is weaker, while at z = the fifth force becomes un- 
suppressed and the density profiles are higher than 
N5 results again for the same reason as above. For 
the most massive halos, the fifth force is strongly 
suppressed even at present, and the halo density 
profiles are lower than in N5 for both z = 1 and 
z = 0. 

Therefore we see the following general evolution pattern 
when comparing the halo density profiles from F and N 



simulations: at early times the fifth force is strongly sup- 
pressed in F simulations, and so the central potential is 
weaker there, producing lower density profiles; then the 
fifth force becomes unsuppressed and the central poten- 
tial becomes as strong as in N simulations, while the par- 
ticles are still moving relatively slowly, producing higher 
density profiles. This happens at z = 1 in F4 simula- 
tions and z = in F5 simulations. Finally the particles 
move fast due to the fifth force, and the density pro- 
files approach the non- Chameleon results as is seen in F4 
simulations at z = 0. The transition happens earlier for 
models with larger \fm\ and for smaller halos (in both 
cases the fifth force becomes unsuppressed earlier). What 
we see in Figs. 12 and 13 are just different stages of the 
above evolution pattern. 

This picture is consistent with the behaviour of the 
power spectrum shown in Fig 6. Once the Chameleon 
mechanism fails and the density profile becomes higher 
than those in N simulations, the power spectrum in F 
simulations tend to be higher on smaller scales (at z = 1 
in F4 simulations and at z = in F5 simulations). Then 
it approaches non-Chameleon results later (at z = in 
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F4 simulations). We expect to see the same evolution 
pattern in F5 and F6 simulations if we would continue 
our simulations in the future. 

In Figs. 12 and 13 we have also shown the efficiency T of 
the Chameleon effect. For the N4 and F4 simulations, we 



see that T = 1 almost perfectly for both z = 1 and z = 0, 
indicating that the Chameleon mechanism does not work 
for \/ro\ = 10 -4 . Note that the agreement between N4 
and F4 serves as an independent test of our code, as the 
treatments for these simulations are very different. For 
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the N5 and F5 simulations, we see that T is close to one 
for the small and medium-sized halos at both redshifts, 
while it is significantly less than 1 for the massive halos, 
especially at z = 1, which are all as expected because the 
Chameleon effect is stronger in high density regions and 
at early times. T is not perfectly equal to 1 even at late 
times and in small halos, because there is a suppression 



of the fifth force and \S/r\ even in the linearised treat- 
ment above the Compton wavelength. In the N6 and F6 
simulations, we see that T <C 1 for all halos and at all 
times, showing a strong Chameleon effect. 

Interestingly, in the F simulations T increases with the 
distance from the halo centre, while in N simulations the 
trend is just opposite (this is clearest in the F6/N6 sim- 
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ulations at z = 0). The former is because the Chameleon 
effect gets weaker as one moves from the halo centre 
(highest density region) and so T tends to 1; the latter is 
because the value of S/r is only affected by particles ly- 
ing within the Compton wavelength from the considered 
position and as one moves from the halo centre more and 
more particles in the central region of the halo become 
unable to affect while this does not happen to $ as 
gravity is a long-range force. 

In summary, we find that the halo density profiles in 
the f(R) gravity model show complicate but interest- 
ing features, which could be understood qualitatively. We 
would like to study these in more details in future works, 
perhaps with the aid of even-higher-resolution simula- 
tions. 



V. CONCLUSION 

Modified gravity scenario, especially the f(R) gravity, 
attracts more and more attention as an alternative to 
dark energy to explain the origin of the cosmic acceler- 
ation at low redshifts. Much efforts have been made to 
investigate the phenomenology of f(R) theories on linear 
scales [6, 10, 22, 23, 34, 35]. Unfortunately, the current 
observational data on linear scales, including weak lens- 
ing, Integrated Sachs- Wolfe (ISW) and cosmic microwave 
background (CMB), etc., can only weakly constrain f(R) 
gravity [36-38]. And even futuristic observations on lin- 
ear scales can hardly prove, or falsify the f(R) scenario 
[35]. 

On nonlinear scales, f(R) gravity can be much better 
tested by cosmological observations [38, 39] because the 
scale dependent enhancement of the growth rate makes 
deviations from GR larger and larger on smaller scales. 
However, we should take into account the mechanism to 
recover GR that is necessary to evade the strong con- 
straints in the solar system. In f(R) gravity models, this 
recovery of GR is accomplished by the Chameleon mech- 
anism. In order to realise the Chameleon mechanism, the 
scalar mode in f(R) gravity should satisfy a highly non- 
linear evolution equation. In this work, we implement 
the Newton-Gauss-Seidel relaxation solver into the MLAPM 
code, an adaptive particle mesh simulation code, and ob- 
tain high resolution simulation results (up to k ~ 20 
h/Mpc for the matter power spectrum). We indepen- 
dently confirmed the results, including the power spec- 
trum and halo statistics, presented in Refs. [15] and [16], 
and extended the resolution by a factor of 7. 

Deviations from GR in the power spectrum are highly 
suppressed at early times when the Chameleon mecha- 
nism works very well. Later the Chameleon starts to fail 
once the background Jr field becomes comparable to its 
fluctuations, 5/r. Then the power spectrum approaches 
the one in the case without the Chameleon mechanism. 
This transition happens at a higher redshift for a larger 
\/ro\ simply because the background field is larger for 
larger |/ho|- We found an interesting behavior that the 



power spectrum in full Chameleon simulations has higher 
amplitude than the one in non- Chameleon simulations 
during this transition. We observed a similar behavior in 
the density profile. The density profile in full simulations 
becomes higher than the one in non- Chameleon simula- 
tions on small scales once the Chameleon has started to 
fail. Qualitatively, this is due to the fact that, once the 
Chameleon mechanism fails, the fifth force becomes un- 
suppressed and the central potential of halos becomes as 
strong as in the non- Chameleon simulations but the par- 
ticles are moving still relatively slowly. Then dark matter 
particles are temporarily trapped at the central region of 
halos. Later particles velocities catch up and the density 
profile and the power spectrum approach those in non- 
Chameleon simulations. In order to confirm this picture, 
we would need higher resolution simulations and study 
the behavior of dark matter particles in halos carefully. 
We can also make the connection between the halo den- 
sity profile and the power spectrum using the halo model 
approach. We leave the study of these issues in a separate 
paper. 



We measured the profile of the scalaron field inside 
halos and studied the efficiency of the Chameleon mech- 
anism As expected, we found that the Chameleon mech- 
anism works better for heavier halos because the den- 
sity is high in these halos. The halo mass function also 
shows the same tendency - the number of heavier ha- 
los approaches that in GR simulations if the Chameleon 
mechanism works as in the \/ro\ = 10 -6 case. For the 
\Iro\ — 10 -4 case, the Chameleon no longer works at 
z — 0, which results in more and more heavier halos 
compared with GR. 



In conclusion, we studied the effect of the Chameleon 
mechanisms in details in our high resolution simulations. 
It was found that once the Chameleon mechanism stated 
to fail, the power spectrum and halo properties showed 
very interesting behaviours before they approach those in 
non- Chameleon simulations. For the power spectrum, we 
showed that this transition can be described by extend- 
ing the Post-Parametrised-Friedmann (PPF) formalism 
that interpolates between the power spectrum in non- 
Chameleon models and GR models. This kind of fitting 
formulae will be useful when we confront our predictions 
to observations as it is far easier to run non- Chameleon 
simulations with different cosmological parameters. Our 
results indicate that \/ro\ = 10 -5 is the critical case 
where the Chameleon mechanism fails to work today. It 
becomes significantly harder to detect deviations from 
GR once \/ro\ becomes smaller than 10 -5 . In this case, 
we need to look into the places where the Chameleon 
still fails to work, i.e. smaller halos. We leave a study 
of observational constraints on f(R) gravity using non- 
linear clustering such as weak lensing and cluster and 
voids abundance in a forthcoming paper. 
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The full equation for u, Eq. (27), contains the quantity 
V 2 e n = V • (e u \7u). To discretise it, let us define b = e M , 
and assumed that the discretisation is performed on a 
grid with grid spacing h. We shall require second order 
precision which is the same as the default Poisson solver 
in MLAPM, and then \7u in one dimension can be written 
as 
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Appendix A: Discretisation of Equations 



where a subscript j means that the quantity is evaluated 
on the j-th point. The generalisation to three dimensional 
case is straightforward. 

The factor b in V • (bVu) makes this a standard variable 
The last thing we are going to do before implementing coefficient prob i em . We need also discretise b, and do it 



the f(R) equations into the A-body code is to discretise 
them, such that they fit to the philosophy of the numeri- 
cal computations - solving the equations on meshes with 
finite grid sizes. This appendix displays the discretised 
equations, and again we discuss the Chameleon and non- 
Chameleon cases separately. We shall only write down 
the discrete equations for u or as those for <£ c as 
simple. 



in this way (again for one dimension) [24]: 
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Then the discrete version of Eq. (27) is 
L h (ui^ k ) = 0, 

in which 
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Then the Newton- Gauss- Seidel iteration says that we can u fj%' usm g our knowledge about the old (and less acu- 
obtain a new (and usually more accurate) solution of u. 
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rate) solution uf^ k as The old solution will be replaced with the new one once 

the latter is ready, using a red-black Gauss- Seidel sweep- 
L h (uf-jj^j m g scheme. Note that 

U iJ% = U i,3,k 7 • (A6) 
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In principle, if we start from some high redshift, then the 
initial guess of 1^,*; could be chosen as the background 
value because we expect that any perturbations should be 
small then. For subsequent time steps we can use either 
the solution at the last time step or some analytical ap- 
proximated solution as the initial guess (for Chameleon 
simulations we find the latter more convenient while for 
non-Chameleon simulations it is just the opposite). 
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2. Non- Chameleon simulations 

The discretisation of the equation Eq. (31) for the non- 
Chameleon simulations is much easier, and similarly to 
the above we have 



L h (v iijtk ) = 0, 
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